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ABSTRACT 

Standard models of accretion discs study the transport of mass on a viscous timescale 
but do not consider the transport of magnetic flux. The evolution of a large-scale 
poloidal magnetic field is however an important problem because of its role in the 
launching of jets and winds and in determining the intensity of turbulence. As a 
consequence, the transport of poloidal magnetic flux should be considered on an equal 
basis to the transport of mass. In this paper, we develop a formalism to study such a 
transport of mass and magnetic flux in a thin accretion disc. The governing equations 
are derived by performing an asymptotic expansion in the limit of a thin disc, in the 
regime where the magnetic field is dominated by its vertical component. Turbulent 
viscosity and resistivity are included, with an arbitrary vertical profile that can be 
adjusted to mimic the vertical structure of the turbulence. At a given radius and time, 
the rates of transport of mass and magnetic flux are determined by a one-dimensional 
problem in the vertical direction, in which the radial gradients of various quantities 
appear as source terms. We solve this problem to obtain the transport rates and the 
vertical structure of the disc. The present paper is then restricted to the idealised case 
of uniform diffusion coefficients, while a companion paper will study more realistic 
vertical profiles of these coefficients. We show the advection of weak magnetic fields 
to be significantly faster than the advection of mass, contrary to what a crude vertical 
averaging might suggest. This results from the larger radial velocities away from the 
mid-plane, which barely affect the mass accretion owing to the low density in these 
regions but do affect the advection of magnetic flux. Possible consequences of this 
larger accretion velocity include a potentially interesting time-dependence with the 
magnetic flux distribution evolving faster than the mass distribution. If the disc is not 
too thin, this fast advection may also partially solve the long-standing problem of too 
efficient diffusion of an inclined magnetic field. 

Key words: accretion, accretion discs - magnetic fields - MHD - ISM: jets and 
outflows - galaxies: jets. 



1 INTRODUCTION 

The presence of a large-scale magnetic field in an accretion 
disc has several potentially important consequences. It may 
play an essential role in the ac celeration and collimation 
of jets and winds from the disc l|Blandford fc Pavnelll982h 
and i n harnessing the rotation al energy of a central black 
hole (fBlandford & Znaick 1977|). It also strongly affects the 
intensity of magnetohydrodynamic (MHD) turbulence and 
the resulting transport o f angular momentum within the disc 
l|Balbus fc Hawlevtll998[ ). 

The origin (and existence) of such a large-scale mag- 
netic field is, however, debated from a theoretical perspec- 
tive. There are in principle two possibilities: either the field 
is created in situ by a dynamo process, or it is brought in by 



the gas that is being accreted. It is still unclear whether a 
dynamo process can create a significant magnetic field that 
would be coherent over a scale comparable to the radius. 
Indeed the MHD turbulence, presumed to be at the origin 
of the dynamo, most likely has a radial correlation length of 
the order of the vertical thickness of the disc, which is usu- 
ally m uch smaller the radius (see the discussion in ISpruitl 
l2010t ). 

In the second possibility, an initially weak field present 
in the gas supplied to the disc could be strongly amplified as 
it is transported inwards and the magnetic flux accumulates 
in the central part of the disc. The transport of magnetic 
flux is due to two processes: the inward advection by the 
accretion flow, and the diffusion due to a turbulent resistiv- 
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ity which tends to counteract the advection and therefore 
transport the magnetic flux outwards. The scenario there- 
fore requires that the field diffuse outwards at a slower rate 
than it is advected for the total transport to be directed 
inwards (at least initially). Which of the advection or the 
diffusion is faster is not obvious a priori, as both are due 
to the same phenomenon: the turbulence which is causing 
an effective viscosity (enabling advection) and an effective 
resistivity (enabling the magnetic field to diffuse). The first 
theoretical studies of the evolution of the magnetic field in 
the presence of advection and diffusion have found that the 
diffusion is much faster than the advection, if t he disc is thin 
and the field lines bend sign ificantly across it (|Lubow et al.l 
11994 iHevvaerts et "afl ll99d ) . This can be understood by the 
following argument. If the angular momentum transport is 
due to a turbulent viscosity v, the advection speed is ap- 
proximately 
3 v 
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The diffusion of the magnetic flux is predominantly due 
to the bending of the magnetic field lines across the disc, 
which is associated with an azimuthal electric current J$ ~ 
B r /{fM)H), where H is the scale-height of the disc. The ac- 
tion of the turbulent magnetic diffusivity r\ on this current 
induces a diffusion speed 

(2) 

Equating the advection and diffusion speeds, one finds the 
maximum inclination of the magnetic field lines that can be 
induced by the advection: 

B r 
B~, 

where V = vjr\ is the turbulent magnetic Prandtl number, 
which is usually expected to be of order unity. The bending 
is thus very small for a thin disc, unless the magnetic Prandtl 
number is unexpectedly large (V ~ r/H). 

This result is rather disappointing as it is problematic 
for models of jets and winds. Indeed, the acceleration of a 
jet or wind by the magneto-centrifugal mechanism requires 
the field lines to bend by more than 30 de grees from the 
vertical direction (|Blandford fc Payne! [l982). Furthermore, 
outward bending of the field lines is a natural consequence of 
the accumulation o f magnetic flu x in th e central part of the 
disc. The result of iLubow et al.l (|l994f ) thus suggests that 
the inward advection of magnetic flux cannot significantly 
amplify the magnetic field, and raises doubts on the very 
existence of a significant large-scale magnetic field. 

Various potential solutions to this magnetic diffusion 
problem have been considered: 

• A high effective magnetic Prandtl number associated 
with the turbulence c ould solve the problem. However, t he- 
oretical expectations l|Parker|[l97ll ; iPouquet et ailll976l ) as 
well as shearing-box simulations of MHD turbulence sug- 
gest that the turbulent magnetic Prandtl number is of or- 
der unity llLesur fc Longarettll2009l ;lG uan fc Gammidl2009l : 
iFromang fc Stone! 12009^ Note however that the main con- 
tribution to the diffusion (the vertical diffusion of a radial 
field) has not been measured so far. It is also possible that 
the effects of turbulence are more complicated than is de- 
scribed by effective diffusion coefficients. 



• B reaking of the axial symmetry: ISpruit fc Uzdenskvl 
(2005) proposed a scenario allowing a fast advection of the 
magnetic flux while preventing its diffusion. In this scenario 
a weak magnetic field concentrates into strongly magnetised 
small patches, owing to the tenden cy of turbulent flows to 
expel magnetic flux fe.g. lWeissll966l ). The diffusion is slowed 
down because the strong magnetic field prevents turbulence 
inside the patch. On the other hand, rapid advection occurs 
as the patch loses angular momentum through a magneto- 
centrifugal wind. This interesting idea needs however much 
more work to be confirmed. 

• Vertical structure: the analysis by ILubow" et al.l (|1994T ) 
uses a vertical averaging of the disc, which implicitly as- 
sumes that the magneti c flux is advected at the same speed 
as the mass. However, I Ogilvie fc Livid (120011 ) have ques- 
tioned the validity of this averaging procedure and sug- 
gested that the vertical dependence of the radial veloc- 
ity, the density and the resistivity could have an influ- 
ence. One obj ective of this article is to inve s tigate furthe r 
this question. Bisnovatvi-Kogan fc Lovelace d2007l. 20121 ') , 



iRothstein fc Lovelace! (|2008l ). and lLovelace et alj(|2009h also 
considered the vertical structure of the disc, studying the 
effect of a non-turbulent layer at the surface. They sug- 
gested that most of the bending of the magnetic field lines 
would take place in the highly conducting non-turbulent 
layer. The low resistivity of the layer would thus drasti- 
cally reduce the diffusion of the magnetic field. A proper 
self-consistent calculation of thi s effect is still m i ssing, how- 
ever, because the solutions o f Lovelace et al.l (|2009h and 
iBisnovatvi-Kogan fc Lovelace! (|2012T l stop at the base of the 
conducting surface layer. The formalism developed in this 
article will be used in a companion paper to test this sce- 



The main purpose of this article and the companion one 
(Guilet & Ogilvie in preparation, hereafter called paper II) 
is to clarify how the vertical structure of the disc affects the 
transport of magnetic flux. We are particularly interested in 
the poloidal magnetic flux threading the disc, as this quan- 
tity satisfies a conservation equation and can be modified 
only by advective or diffusive transport, possibly enhanced 
by turbulence. Since the poloidal magnetic field strongly af- 
fects the intensity of turbulence in an accretion disc and is 
an essential component in magneto-centrifugal jet launch- 
ing, the transport of poloidal flux needs to be considered on 
an equal basis to the transport of mass in an accretion disc. 

Although future work may involve direct numerical sim- 
ulations of turbulent discs, in our present analysis the tur- 
bulence is simply modelled by an effective viscosity and an 
effective resistivity. The formalism described in Sections [2] 
and Section [3] allows any vertical profile of these coefficients, 
which will be used in paper II to study the effect of the 
variation of the turbulent intensity with distance from the 
mid-plane. The rest of this article however focuses on the 
idealised case of uniform diffusion coefficients. Under these 
assumptions and in the limit of a thin accretion disc, we self- 
consistently solve for the vertical profile of the mean velocity 
and magnetic field. We focus in particular on the regime in 
which the magnetic field is dominated by its vertical com- 
ponent. This is done partly to simplify the analysis and take 
advantage of a linear system of equations, and partly in or- 
der to avoid the complications of magneto-centrifugal jet 
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launching from th e disc. A sepa rate paper will address the 
launching process (|Ogilvidl20ld) . 

The plan of this paper is as follows. In Section[2] we use 
the basic equations of MHD, together with an asymptotic 
expansion in the limit of a thin disc, to derive the equations 
describing the evolution of the mass and magnetic flux dis- 
tributions on a viscous or resistive timescale. In Section O 
we focus on the quasi-local problem describing the vertical 
structure of the thin disc and enabling the calculation of the 
transport rates. This problem is solved for the special case 
of uniform diffusion coefficients in Section]!] The results are 
discussed and summarised in Section [S] 



2 AN ASYMPTOTIC EXPANSION WITHIN A 
THIN DISC 

We start from the equations of MHD, in the form 

dp 



at 



+ V ■ {pv) = 0, 



dv 



+ v ■ Vd = -pV$ - Vp 



1 



OB 

dt 



H (V x B) x B + V ■ T, 

Po 



= V x (v x B - r]V x B), 



V-B = 0, 



pv 



Vu + (Vu) 1 - -(V-w)I 



(4) 

(5) 
(6) 
(7) 

(8) 



where p is the density, v the velocity, $ the gravitational 
potential, p the pressure, B the magnetic field, T the viscous 
stress, r\ the magnetic diffusivity, v the kinematic viscosity 
and I the unit tensor of second rank. (A bulk viscosity could 
be included, but would have no effect on this problem.) We 
neglect self-gravity and do not consider a thermal energy 
equation at this stage. 

We consider an axisymmetric solution of these equa- 
tions in cylindrical polar coordinates (r,(f>,z). The poloidal 
part of the magnetic field can be described in the usual way 
by a poloidal magnetic flux function ip(r,z,t). Thus 



plane z = 0. In the absence of magnetic fields, there is a 
well understood ordering scheme for thin discs. If the char- 
acteristic angular semithickness of the disc is H/r — O(e), 
where < e < 1 is a small dimensionless parameter, and 
the Shakura-Sunyaev viscosity parameter is a, then the ra- 
tio of the sound speed to the orbital velocity is O(e), and the 
ratio of the radial velocity to the orbital velocity is O(o.e 2 ). 
The fractional deviation of the azimuthal velocity from the 
Keplerian value is 0(e 2 ), and so on. Furthermore, the disc 
evolves on the viscous timescale, which is longer than the or- 
bital timescale by 0(a^ 1 e~ 2 ). Although this is rarely done, 
the equations governing thin accretion discs can be obtained 
by a formal asymptotic expansion of the basic equations of 
fluid dynamics. 

For magnetised discs there is more than one possible or- 
dering scheme, depen ding on the st rength (and orientation) 
of the magnetic field (|Ogilvidl 19971 ). We consider here a sit- 
uation in which the disc has comparable values of v and 77, 
with (formally) a = O(l). The magnetic field is dominated 
by the vertical component, which has a magnetic pressure 
comparable to the gas pressure. This assumption allows us 
to avoid a consideration of magnetocentrifugal jet launching. 

We resolve the internal structure of the thin disc and its 
slow evolution in time through the introduction of a rescaled 
vertical coordinate 



c = e - 1 * 

and a time variable 



t = at. 



(11) 



(12) 



We then propose the following expansion of the fluid vari- 
ables. (The scheme for indexing the variables is debatable, 
but the choice made here is convenient for the present pur- 
poses.) 



P = Po(r, C, t) + <?pi(j, C, t) + 0(e 4 ), 



P = <? [po(r, C, t) + € 2 p 2 (r, C, r) + 0(e 4 )] 



$ = $ (r) + e 2 $ 2 (r)iC 2 + O(e 4 ) 



v r = e 2 v r2 (r, C, t) + e 4 v r4 ,(r, (, r) + 0(e 6 ), 



(13) 



(14) 



(15) 



(16) 



-e<f, x Vi[> + B$ e^. 



(9) 



This representation satisfies the constraint V ■ B = auto- 
matically. Poloidal magnetic field lines correspond to curves 
ip — constant in the (r, z) plane, and (up to an additive con- 
stant) 2nif) is the magnetic flux contained within a circular 
loop at position (r,z). The poloidal part of the induction 
equation can then be integrated to give 



dt 



+ v ■ Vtf> 



-T-7 

rjr V • 



(10) 



which shows that the poloidal magnetic flux satisfies a con- 
servation law and is affected by advection and diffusion. 

We are interested in solving these equations both in- 
side and outside a thin accretion disc that lies close to the 



W0 = rOo(r) + e 2 v<t,2(r, £, r) + e^v^r, £, r) + 0(e 6 ), (17) 



v z = e [e v z2 (r, <, r) + e v z4 .(r, f, r) + 0(e 



^ = e[Mr,r) + e 2 Mr,C,r) + 0(e 4 )] , 



= e [eB^i (r, <, r) + e 3 B 03 (r, C, r) + 0(e 5 



(18) 



(19) 



B r = e [eB rl (r,C,r) + e 3 B r3 (r,C,r) + 0(e 5 )] , (20) 



(21) 
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B z = e[B z0 (r,T) + e 2 B z2 (r,Cr) + 0(e 4 )] , 



(22) 



= B z0 d ( v r2 + d ( [r)o(d ( B r i - d r B z0 ) 



(30) 



u = e 2 [u (r,C,r) + O(e 2 )], (23) 



7 ? = e 2 [7 ?0 (r,C,r) + O(e 2 )]. (24) 

We make the following observations. The overall scaling of 
the density is irrelevant provided that self-gravity is ne- 
glected; however, the relative scaling of the pressure and 
density is significant and implies that the sound speed is 
O(e) (in comparison to the orbital velocity). The expansion 
of $ is the Taylor expansion of a smooth, symmetric ex- 
ternal potential about the midplane z = 0; in the case of 
a point-mass potential $ = —GM(r 2 + z 2 )~ 1 '" 2 , we have 
$o = -GM/r and $ 2 = GM/r 3 . The scaling of B is such 
that, as mentioned above, the magnetic field is dominated 
by the vertical component, which has a magnetic pressure 
comparable to the gas pressure. The expansion of ip is of 
the form required to produce this. The scalings of v and rj 
correspond to having (formally) a — O(l). 

Substituting these expansions into the basic equations 
and comparing the coefficients of powers of e, we obtain a 
succession of equations. The radial component of the equa- 
tion of motion at leading order yields 

-porfl 2 , = -p d r ®o, (25) 

which implies rfi 2 , = d r &o, i.e. that Qo(r) is the angular 
velocity of a circular particle orbit of radius r in the mid- 
plane. In the case of a point-mass potential, we obtain the 
Keplerian value Q = (GM/r 3 ) 1/2 . 

The vertical component of the equation of motion at 
leading order yields 



= Brir9j.no + B z0 9 c t^2 + d^rjod^B^) 



(31) 



= -po*2C - %J , 



(26) 



which is the usual condition of hydrostatic equilibrium in 
which the vertical pressure gradient balances the vertical 
gravitational force. Under the present scalings, the Lorentz 
force does not affect this balance at leading order. The so- 
lution of this equation depends on additional assumptions 
regarding the thermal physics of the disc. In the simplest 
situation of a vertically isothermal disc in which po = c^oPo, 
where the isothermal sound speed c s o(r, r) does not depend 
on £, we obtain the familiar solution 

*2 



PO = 



(2n) 1 / 2 Ho 



exp 



2H 2 



(27) 



where Ho(r,r) = c s o/&/ 2 is the isothermal scaleheight and 
So (r, r) is the surface density. 

Next, the horizontal components of the equation of mo- 
tion and of the induction equation yield the four equations 

32 



-2p o n o V02 = -pod r <&2^C 
B z o 



H d c B rl + d ( (p vod ( v r 2), 

Po 



Po«r2-9 r (r 2 n ) 

r po 

+d i (povod(V 4>2 ), 



" z0 9 c B^,i + ^d r {povor 3 d r £l Q ) 



(28) 



(29) 



These can be regarded as four linear equations for the un- 
knowns v r 2, i^2, B r i and B^i. They are linear because of 
the assumption of small deviations from orbital motion and 
a vertical magnetic field, implicit in the asymptotic expan- 
sion. Although these four quantities depend on r, £ and r, 
only derivatives with respect to £ appear in these equations 
and they can therefore be regarded an eighth-order system 
of ordinary differential equations (ODEs) in £ at each r and 
r. This situation arises because the vertical dependence is 
assumed to be more rapid than the radial dependence in a 
thin disc, and the time-dependence is assumed to be slow. 
The equations are inhomogeneous and are driven by various 
source terms: the vertical dependence of the radial gravi- 
tational force, the radial gradient of total pressure, the az- 
imuthal viscous force resulting from the orbital motion, etc. 

The expected symmetry of the solutions is that v r 2 and 
v,j,2 are even functions of £ while B r i and B^i are odd. The 
behaviour of the solutions at large |£| can be illustrated by 
considering the case of a vertically isothermal disc in which 
uo and r/o are bounded as £ — > oo (e.g. if they are indepen- 
dent of £). In this case the solution is of the form 



Vr2 



Ci + £i(C), 



V4.2 ~ C 2 ±< 2 ± C 3 C + C* 4 + E 2 (0, 



Bri ~C 5 {±C 6 + E 3 (<;), 



B^~±C 7 + E 4 (() 



(32) 



(33) 



(34) 



(35) 



as £ — ► ±co, where the d are constants and the Ei are 
(like the density and pressure) exponentially small; here the 
parametric dependence of the solution on r and r has been 
suppressed. In order to satisfy the governing equations, the 
values of the constants are constrained by 



C2B z0 = ~Cr>rdM), 



CzB z o 



Cr, = d r B z o; 



(36) 



(37) 



(38) 



while C@ and Cr have the nature of input parameters, Ci 
and C4 have the nature of output parameters. The C5 term 
is required in order to produce a force-free field in the low- 
density exterior at large £|. (In fact it is a current-free or 
potential field, because an axisymmetric poloidal force-free 
field must be current-free.) The Ce term represents, in some 
sense, the inclination of the poloidal field at the upper sur- 
face of the disc, and we call G\ = B rB accordingly. More 
precisely, it is the radial component of the magnetic field ob- 
tained by extrapolating the parabolic form of the field lines 
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in the force-free region above the disc down to the midplane 
£ = 0. The CV term represents the azimuthal component of 
the magnetic field at the upper surface and we call it B^. 
If we are considering discs without outflows (because of the 
condition |B rs | <C \B Z \) then we should set £?0 S = in order 
that the magnetic stress B^Bz/fio vanish at large £j. How- 
ever, we retain this term as a convenient way to model the 
effect of angular momentum removal by an outflow with- 
out the complications of solving for the acceleration of the 
flow. Angular momentum is removed by specifying a nega- 
tive value for B t j >s B z . 

In terms of the flux function ip, we also have the rela- 
tions 

B z0 = -dri>o, (39) 



Brl = 



1 



(40) 



Neither tpo nor B z o depends on so B z o can be regarded 
as an input parameter in the above system of ODEs. 

Furthermore, the integrated form of the poloidal part 
of the induction equation (equation [10]) yields 



9-rVo + v r2 d r ipo = Vo 



rd r (J^dripo^j + d%ip2 



(41) 



Note that — (l/r)9^ of equation (|41|) produces equation (|30l) . 
Once we have solved the system of ODEs, equation (|41|) can 
be evaluated at any value of ( (or indeed averaged in any 
convenient way with respect to £) to discover the value of 
d T ipo- This gives the rate of poloidal magnetic flux transport 
at the given values of r and r , which we can associate with an 
effective transport velocity Vtj, through d T tpo + v^dr^o = 0, 
which implies d T B z o + (l/r)d r (rv^,B z o) = 0. In fact, if we 
evaluate equation (|41|) in the limit £ — > oo, we find dripo = 
—CirB z o and therefore v^, = G\. Since the electric current 
vanishes in this limit, the radial velocity corresponds to the 
speed at which the magnetic field is transported. 

In addition, the equation of mass conservation yields 



d r Qo, d^Oo, d r Cso, <9r-£o, d r vo, 9 r r]o and d r B z o- The output 
parameters of greatest interest are v m and v^,. The problem 
can be simplified, and made dimensionless, by assuming a 
Keplerian disc (Oo oc r~ 3 ^ 2 and $2 = Oq) and alpha pre- 
scriptions of the form uo = ac 2 /Qo, Vo = vq/V with a and 
V being constants. 

In this analysis we have simplified the problem by as- 
suming that the magnetic field is dominated by its vertical 
component within the disc. This field needs to be matched 
to a force-free (in fact, current-free) field outside t he disc, 
which involves the solution of a global problem (e.g. [OEilvic 
[1997]). The nature of this problem is that the poloidal flux 
distribution on the disc determines the inclination of the 
poloidal field at each radial location, and in general this 
angle is not small. An alternative ordering scheme is pos- 
sible, in which all three components of the magnetic field 
are O(e). In such a scheme, the inclination of the poloidal 
field could be large and magneto-centrif ugal jet launching 
would be possible l|Ogilvie fc Livid feoOll ). The pressure of 
the horizontal magnetic field would also contribute to the 
vertical force balance and the problem would become more 
nonlinear. (In fact, it is easy to add this effect by hand to 
the present system of equations, but we will not do so here.) 
However, the transport of magnetic flux by diffusion would 
be formally faster by a factor 0(e _1 ) than that by advec- 
tion, if V = O (l); this ineq u ality is merely an expression 
of the result of iLubow et alj l| 19941 '). Therefore no ordering 
scheme is fully satisfactory. We adopt the present scheme 
because of the transparency of its derivation, the linearity 
of the resulting equations, and because it avoids the compli- 
cations of jet launching, while missing only the effect of the 
magnetic compression of the disc. 



d T pa + -dr{rpov r 2) + d c (pov z2 ) = 0. 



(42) 



We define the surface density at leading order, 

/oo 
po(r,C,r)dC. (43) 
- OC 



Then (assuming no vertical mass loss from the disc) 



(44) 



9r£o + ~dr \T J pOVr2 d( j = 0. 

Let 

mo(r,r) = / S (r', r) 2-kt' dr' 
Jo 

be the mass contained within radius r at time r; then 9 T mo + 
v m d r mo = 0, where 



(45) 



«m = ^T / P0Vr2 d( 
^0 . 



(46) 



is the effective transport velocity for mass, which can be 
determined from the solution of the system of ODEs. 

In summary, then, the input parameters at given val- 
ues of r and r are O , $2, c a0 , S , v , r/ , B z0 , B rs , B^, 



3 A QUASI-LOCAL PROBLEM 



Equations (|28|) (|31[) define a quasi-local problem, as only 
vertical derivatives of the unknowns appear, while radial 
derivatives appear only as source terms. The solution of 
this local problem provides the radial transport velocities 
of mass and poloidal magnetic flux (at a given time and 
radius), which are necessary to determine the global evolu- 
tion of the disc on a viscous/resistive timescale. The local 
problem on the other hand is effectively stationary (no time 
derivatives appear), because stationarity is assumed on the 
much shorter dynamical timescale. In this article and pa- 
per II, we solve this local problem and leave for future work 
the study of the global evolution of a disc. In this section 
we rewrite in a dimensionless form the equations governing 
the local problem. We discuss the dependence on the differ- 
ent parameters, and describe the numerical method used to 
solve the equations. 
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3.1 Differential system 

We choose as independent variables po, v r 2, w^i, B r i and 
B^i. These variables are non-dimensionalised by defining 



ters: 



b,j, 



PoH 

E ' 

r Vr2 

H c s ' 
r v<t,2 
H c s ' 
r B r i 
H B z ' 
r B^i 
H B z ' 



(47) 
(48) 
(49) 
(50) 
(51) 



where we have dropped the subscript on quantities such as 
H , Cs, B z and E. A dimensionless vertical spatial coordinate 
is also defined as: 



C = z/H. 



(52) 



Note that this variable is different from £ used in Section [5] 
which was a rescaled but dimensional variable. In the rest 
of the paper when we use the symbol £, we refer to the 
dimensionless variable. 

We assume a point-mass potential and therefore circu- 
lar Keplerian orbital motion at leading order. We use the 
standard alpha prescription for the viscosity, v = ac^/Q. = 
ctCsH, allowing for a vertical (but not radial) dependence of 
the a parameter through a = aog(£). The resistivity r\ is 
then related to the viscosity through the magnetic Prandtl 
number V = v/r), which can also have a vertical dependence. 
An isothermal equation of state is assumed for simplicity. 
Equation (|27|) gives the corresponding density profile, which 
can be written in dimensionless form as 



1 



-C 2 /2 



(53) 



The differential system given by equations ([28}-([3T} is 
rewritten as 

1 13 

- -de (padtUr) - 2us - -—d c b r — - + D H - A,s 
p p p 2 



\dc [pad^u^,) + \u r - -J^debd, 
p 2 fj p 



(54) 



3 

2° 



+D H 1 + 



din a 
dln< 



(55) 



(56) 



D B 



v 2 

Po ^e 3 
Bl H 
d\nH 
d\nr ' 

2D H - 
d\nB z 



3 <91nS 
2 + 91nr ' 



(58) 
(59) 
(60) 
(61) 



<91nr 

Here /3o corresponds roughly to the midplane value of the 
plasma ft parameter (the ratio of the thermal pressure to 
the magnetic pressure); more precisely, the two are re- 
lated by /3(£ = 0) = \/2/-k Po- The parameter D„£ equals 
<91n(^E)/(91nr, given that v = aH 2 Q and a does not de- 
pend on r. 



3.2 Boundary conditions 



From equations (|32|l (|35[) , one can deduce that the following 
quantities vanish exponentially fast as £ — > ±oo: 



pu r 



0. 



pu$ — > 0, 



{D B (±b r 



0. 



(±6 



0. 



(62) 



(63) 



(64) 



(65) 



Let us recall that these conditions are motivated by the ab- 
sence of outflow, owing to the low inclination of the field 
lines. As a consequence the magnetic field at infinity is force- 
free. The inclusion of a non-vanishing azimuthal component 
of the magnetic field is not self-consistent but allows us, if 
we wish, to mimic the effect of angular momentum removal 
by an outflow. These boundary conditions are homogeneous, 
except for the linear source terms proportional to Db, b rs , 
and bfa in equations (|64[) and (I65|) . 

As the differential equations, the source terms and 
boundary conditions have reflectional symmetry about the 
midplane, the stationary profile we seek shares the same 
symmetry (implying a symmetric horizontal velocity and an 
antisymmetric horizontal magnetic field). As a consequence, 
the following conditions apply at the midplane ( — 0: 



d(U r = 0, 



d^Utf, = 0, 



K = 0, 



(66) 



(67) 



(68) 



-d< (^d c b^j - d<Uf + ^b r = 0, (57) 
where we have defined the following dimensionless parame- 



0. 



(69) 



Reflectional symmetry also implies that the solution need 
be computed only in the half-space £ > 0. 
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3.3 Transport velocities of mass and magnetic flux 

The transport rates of mass and magnetic flux can be ob- 
tained from the solution of the above equations and ex- 
pressed in terms of a transport velocity in the following way. 
Multiplying the azimuthal component (|55[) of the equation 
of motion by 2p, integrating over all £ and applying the 
boundary conditions, we obtain 

/ ( 5« r dC = --(l + 2A, s ) / pad(+ —b^s. (70) 

J-oc, 1 J-oa PO 

In the case of uniform a this gives the transport velocity 

(71) 



3 4 
u m = --a(l + 2D vT .) + —6, 



where 



H c s 



pu r d£ 



(72) 



is the dimensionless mass transport velocity. More generally, 
a should be replaced by oto in equation (|71[) and a factor of 



r 



9(0 dC 



d< 



(73) 



should be included. For a vanishing these results are 
consistent with the familiar expression from the standard 
theory of a Keplerian accretion disc, 

^-^^ (74) 

where the overbar refers to a density-weighted vertical av- 
erage. 

The integrated version of the radial component (|56|) of 
the induction equation is 



— u r + -p(d(b r — Db) = const, 
where 



r ty, 
H c s 



(75) 



(76) 



is the dimensionless magnetic flux transport velocity. It can 
be evaluated using equation (|75p at any convenient hight. 



3.4 Relation to previous works 

lOgilvie fc Livid (|200lft solved for the vertical structure of 
a magnetised accretion disc, allowing for the possibility of 
magneto-centrifugal j et launching (see also lOgilviel 1 19971 ; 
lOgilvie fc Livid 1 19981 ). The equations they solved are sim- 
ilar to ours, but with the following differences. The ther- 
mal structure of the disc was determined by a balance be- 
tween viscous and resistive heating and radiative cooling, 
rather than being assumed to be isothermal. This solution 
was matched to an isothermal atmosphere with a force-free 
magnetic field. The radial gradients of thermal and mag- 
netic pressure were neglected, as was the vertical transport 
of momentum by viscosity. The pressure of the horizontal 
magnetic field was taken into account in the vertical force 

bal ance. 

Lovelace et al. (|2009j) and 

iBisnovatvi-Kogan fc Lovelacel (|2012h also performed a 
similar calculation of the vertical structure of a magnetised 
accretion disc. The differences between their approach and 



ours are the following. The equations in the limit of a thin 
disc were derived more systematically in this paper, and 
describe more p recise ly the vertical structure: contrary to 
I Lovelace et al.l (|2009h we take into account the vertical 
dependence of the pressure gradient and viscous stress. 
We also did not discard the radial gradient of vertical 
magnetic field, whose contribution to the magnetic flux 
transport appears at the same order in the expansion 
as the advection by an accretion flow. Furthermore, we 
solve a more general problem, as we do not assume that 
the disc is stationary on a viscous/resistive timescale but 
only on a dynamical timescale, and we consider a much 
wider range of /3. The vertical boundary conditions also 
differ substantially. We calculate the vertical profiles up 
to the region that is magnetically dominated, where one 
can impose the magnetic field to be th at dictated by 



Lovelace et al.l 
2012ft do not 



the ex terior solution. On the other hand , 
|2009T ) and IBisnovatvi-Kogan fc Lovelacel 
solve the transition between the interior of the disc and the 
magnetically dominated region; instead they apply special 
boundary conditions at the surface of the disc. The validity 
of these boundary conditions may be questioned in light of 
the results of Section [4j since we find a significant jump in b r 
and at the transition to the force-free regime, while their 
boundary conditions preclude such jumps (note however 
that their diffusion coefficients vanish in the surface layer, 
which is not considered in this article but delayed to paper 
II). Another difference is the choice of a: considering the 
transition to the magnetically dominated regime forces us 
to determine a through the marginal stability hypothesis in 
order to avoid unphysical effects (see Section [3.7[) . Finally, 
they use a disc model where the density is independent of 
height inside the disc, while we use an isothermal model 
where the density is a Gaussian function of height. Allowing 
a variation of the density is essential in our result that 
the advection speed of mass and magnetic flux can differ 
dramatically (see Section [4}. 



3.5 Relation to the shearing sheet 

Both our equations and those of lOgilvie fc Liviol (|200ll ) are 
related to the local approximation (shearing sheet or shear- 
ing box), commonly used in the study of accretion discs. 
Starting with the equations of MHD in this approximation 
and assuming a solution that depends only on the vertical 
coordinate z, we would obtain equations equivalent to ours 
except that the source terms would be absent and the pres- 
sure of the horizontal magnetic field would be taken into 
account in the vertical force balance. 



3.6 Relation to the magnetorotational instability 

The four equations (|54[l - (|57p governing the horizontal com- 
ponents of the velocity and magnetic field are closely related 
to those appearin g in the analysis of the magnetorotational 
instability (MRI; IfBalbus fc Hawlevlll998l '). Let us write the 
equations in the symbolic form 



LX 



(77) 



where X = [u r b r b ( j > ] T is a vector of unknowns, L is 
the linear operator that generates the left-hand sides of the 
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equations, and F represents the source terms on the right- 
hand sides. To analyse the MRI of an isothermal disc with a 
vertical magnetic field, we would instead solve the eigenvalue 
problem 

LX = XX, (78) 

with homogeneous boundary conditions, corresponding to 
solutions of the linearized equations proportional to e At and 
with no horizontal spatial dependence. The MRI would cor- 
respond to a solution with Re(A) > 0. The close connection 
between the MRI and the equilibrium of magnetised disc s 
has been noted before |Ogilvidll998l ; IOgilvie fc Lividl200ll i. 

3.7 Marginal stability hypothesis 

The value of a (or ao if a vertical dependence is permit- 
ted) can be estimated by supposing that the MRI is made 
margi nally stable by the eff ect of the viscosity and resis- 
tivity jOgilvie fc Lividl2001i ). This approach has the desir- 
able property of ensuring that the obtained stationary so- 
lution is not unstable (at least to modes with no horizon- 
tal dependence, which are invariably the most dangerous). 
I Ogilvie fc Livid (|200lT ) showed that this assumption avoids 
unphysical effects like multiple bending of the field lines (see 
also Appendix A). 

Stability depends on the strength of the magnetic field 
(through Po) and on the viscosity and resistivity (through 
a and V). As expected, the most difficult mode to stabilise 
is the largest-scale c hannel mode, which is antisymmetric 
about the midplane (|Ogilvie fc Livid 1200 j ). and this mode 
has A = at marginal stability. Q For weak fields (/3 > 10 4 ) 
the marginal mode is localised around ( = ±Cb (Figure [T]), 
where 

Cb = ^(~/3o 2 ) (79) 

is the height at which the magnetic pressure equals the ther- 
mal pressure. 

While physically motivated, this way of determining the 
effective turbulent diffusion coefficients should be taken with 
a grain of salt. Particularly worrying is the tendency of this 
method to produce large values of a. Indeed, in the case 
of uniform diffusion coefficients, a is of order unity over a 
large range of magnetic field strength (Figure[2|). Such a high 
value is most probably not realistic for weak magnetic fields 
{Po 1), although we note that direct numerical simula- 
tions of discs with vertical gravity and a net vertical mag- 
netic flux are computationally difficult and few results are 
available, perhaps precisely because the turbulence is very 
intense when /3o is not very large. The somewhat surpris- 
ingly slow decline of a at large /3o is due to the fact that its 
value is determined mostly by the region where the chan- 
nel mode is localised, which is rather strongly magnetised 
(P ~ 1) even when the midplane is very weakly magnetised. 
Obviously, assuming the same value of a in these regions 
with vastly different magnetisation is unrealistic. The use 

Although this means that the linear operator L is singular and 
possesses an antisymmetric null cigcnfunction, equation 1 1771) still 
has a unique solution that is symmetric about the midplane, be- 
cause the operator is invertible within the symmetric subspace. 



5 



4 r 

3 r 

a : 

2 r 




10° 10 1 10 2 10 3 10 4 10 5 

Figure 2. The value of the a parameter determined by the 
marginal stability hypothesis as a function of the magnetic field 
strength. The calculation is done with uniform diffusion coeffi- 
cients. The full line illustrates the value of a obtained with the 
magnetic Prandtl number mostly used in this paper, V = 1. 
The dotted line shows the a values obtained without viscosity 
(P = 0), and the dashed line the corresponding analytical esti- 
mate for large /3q of Appendix A (equation I149H . 

of a vertical profile for a, as is done in paper II, alleviates 
this issue by reducing the midplane value (q(C = 0) ~ 0.1 is 
more realistic although still probably too high for very weak 
fields), while keeping similar values far from the midplane. 

Fortunately (given the uncertainty on the value of a), 
the balance between the advection and diffusion of the mag- 
netic field is only weakly dependent on a, when the angular 
transport enabling advection is caused by turbulence. In- 
deed, both advection and diffusion processes are (at least 
approximately) proportional to a, through respectively the 
viscosity and the resistivity. 

3.8 Parameters and form of the solution 

The free parameters of the dimensionless problem are the 
magnetic Prandtl number V and the magnetisation param- 
eter Po, as well as the source terms Dh, D v ^, Db, b rs and 
b(f,s- 

An important consequence of the linearity of the equa- 
tions is that the solution depends linearly on the source 
terms, appearing either on the right-hand side of the system 
of equations or as a non-vanishing boundary condition at 
infinity (b rs and 6</, s ). Thus, the general solution is a linear 
combination of the solution vectors corresponding to each 
source term: 

X = X K + XdhDh + X DvT,D v Y, + XdbDb 

+Xf lrs b rs + Xfj^b^s, (80) 

where Xdh is the solution vector corresponding to the 
source term proportional to Dh and so on. Xk is the so- 
lution vector corresponding to the source terms F that are 
independent of Dh, D v t,, Db, b rs , and b<f, s ; these terms con- 
tain the radial derivative of the leading-order angular veloc- 
ity which is assumed to be Keplerian (hence the notation), 
as well as a term describing the vertical dependence of the 
gravitational potential and other geometrical terms coming 
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z/H z/H 



Figure 1. Marginal MRI channel modes for different strength of the vertical magnetic field : /3rj = 10 (red), /3o = 10 4 (blue), and fio = 10 s 
(black). The calculation is done with uniform diffusion coefficients and V = 1. The modes are normalised such that their maximum radial 
velocity equals one. The vertical dotted lines illustrate the height at which the magnetic and thermal pressures coincide, and where 
the marginal channel modes tend to localise. 



from differentiating factors depending on r in the viscous 
term and the radial pressure gradient. Furthermore, we de- 
fine 

Xh y d = Xk + Xdh, (81) 

being the solution corresponding to the hydrodynamic 
source terms with the standard parameters Dh = 1 (rele- 
vant to a disc with a constant aspect ratio H/r) and D v t. = 
(relevant to a steady accretion disc far from the inner bound- 
ary). For these parameters, we thus have 

X — -X^hyd + XdbDb + Xt rs b rB + X^sb^g. (82) 

The linearity of the problem makes the exploration of 
the parameter space and the physical understanding of the 
solutions much easier. Indeed, given the marginal stability 
hypothesis, the solution depends in a nonlinear way only 
on the two parameters j3o and V. For each pair of values of 
these two parameters, one needs to compute the six solution 
vectors Xk, Xdh, XduS, Xdb, X brs and X b<tya , each of 
which can be represented by plotting the profiles of u r , u^, b r 
and b^. The general solution is then just a linear combination 
of these solution vectors with the appropriate coefficients. 

Note that another input to the model is the vertical 
profile of the diffusion coefficients, the shape of which can 
be freely imposed (although its normalisation is determined 
by the marginal stability hypothesis). In Section]!] we study 
in detail the simplest case of a uniform resistivity and vis- 
cosity. The effect of the vertical structure of the diffusion 
coefficients will be considered in paper II. 

3.9 Method of numerical solution 

We solve the problem described in the previous subsections 
with the use of two different methods: the shooting method, 



and a (more successful) spectral method using a decompo- 
sition on a basis of Whittaker cardinal functions, whi ch are 
well suited to problems on an infinite interval (e.g. iBovdl 
l200ll ; lLatter et al.|[20Tol ). In the shooting method, we shoot 
from the midplane by guessing the values of u r , u^,, c\& r , and 
dc^bfj, there (equations 1661 [69]) . The differential system given 
by equations (|54p - (|57[) is then used to evolve this solution 
up to a height £n Um , where the boundary conditions at infin- 
ity are applied. Newton iteration is used to converge to the 
desired solution. £ llum should be chosen large enough to lie in 
the force-free regime, in which case the solution is indepen- 
dent of ("num. We found that choosing ( nnm typically one or 
two scale-heights above C,b was sufficient for this purpose. 
Similarly, the 'grid' of the spectral method should extend 
up to several scale-heights above (b to obtain converged re- 
sults. The decomposition on Whittaker cardinal functions 
implicitly imposes the condition that the variables tend to 
zero exponentially fast at infinity. For this reason, in the nu- 
merical calculation using the spectral method we replace b r 
and btf, by the following variables: 

b r = b r - D B C - Ws tanh(C 3 ), (83) 

fe = bet, - fe 0s tanh(C 3 ), (84) 

which should vanish exponentially fast at infinity accord- 
ing to the boundary conditions stated in equations (1641 1 
(|65[1 . The differential system had to be changed accordingly. 
We obtained very good agreement between the two methods 
(Figure EI . 

We also compared the results (with no source terms 
except b rs ) with one-dimensional time-dependent direct nu- 
merical simulations of a stratified shearing box in which 
the value of B r is imposed at the vertical boundaries. 
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Figure 3. Comparison of ID time-dependent direct numerical simulations of a stratified shearing box (full lines) and solutions by the 
spectral method (dashed line) and the shooting method (dotted line). The black full line is a simulation with B rs = 10~ 3 B Z (it is 
almost indistinguishable from the dashed and dotted lines), while the orange line is a simulation with B rB = 0.5B Z . Other parameters 
are /3q = 10 4 , a = 1.5, V = 1, and all source terms except b rs are set to zero. 



The simulations were performed wit h the code RAMSES 
|Tevssierl |2002| ; iFromang et alj [20061 ). The comparison al- 
lows us to evaluate the limitation due to our assumption 
that \B r \ <C \B Z \ (Figure [3]). Fortunately even when the 
field is significantly inclined (B r = 0.5S Z , just below the 
inclination threshold above which a magnetocentrifugal jet 
is launched), the velocity and magnetic field profiles remain 
very close to those obtained under the assumption of a small 
inclination. This is observed to be true as long as pa is rather 
large. For /3q values of order unity magnetic compression (ne- 
glected here) would be expected to play a significant role if 
the inclination is large. 

The results described in the remainder of this paper 
were obtained with the spectral method, which gave a good 
convergence on a broader parameter space than the shooting 
method. 



4.1 Approximate analytical model 

Before discussing the numerical solutions of our problem, 
we develop an approximate analytical model that is help- 
ful in their interpretation. The magnetic field and velocity 
profiles can be easily understood in the two limiting cases 
of weak (passive) magnetic field (/? > 1), and very strong 
(force-free) magnetic field (/3 -C 1). In the case where the 
magnetic pressure is small compared to the midplane pres- 
sure (Po 3> 1), both regimes appear: close to the midplane 
the field is passive, while at £ — > 00 the field is force- free. The 
transition between the two regimes takes place around (b 
where the magnetic pressure is comparable to the thermal 
pressure (£a is given by equation {7SJ). In the following two 
subsections we describe the two different regimes. In a third 
subsection we construct a simple model that connects the 
two regions to obtain an analytical estimate of the profiles, 
and in particular of the transport velocity of the magnetic 
flux. 



4 THE CASE OF UNIFORM DIFFUSION 
COEFFICIENTS 

The formalism developed in Sections [2] and [3] allows in prin- 
ciple for a vertical dependence of the diffusion coefficients. 
However, in the remainder of this paper we focus on the case 
of uniform diffusion coefficients. This simple, but probably 
unrealistic, model is a first step and has the advantage of 
allowing for some analytical treatment. More general mod- 
els involving non-uniform diffusion coefficients are studied in 
paper II. In the remainder of this section, unless otherwise 
noted, the numerical calculations use our fiducial parameters 
V = 1 and p = 10 4 . 



4-1.1 Passive magnetic field 

In the limit of very large P, the Lorentz force is negligible 
and the velocity profile is unaffected by the magnetic field. 
If the kinematic viscosity is uniform, the purely hydrody- 
namic velocity profiles can be computed analytically and 
are parabolic: 

U r = U r o+Ur2( 2 , (85) 

u,/, = u^o+u^C 2 , (86) 
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with the following coefficients: 



4-1.3 Two-zone model 



U r 2 



U<j>Q 



1 +4a 2 



(3-5-Dh) 



(87) 



n 



l + 4a- 



- 2 \D H 



(3 - 5D H ) 



1 +4q ; 



(3-5Dh), 



^4>2 = 



2 V 2 



+ 



1 + 4a 2 



(3 - 5D H ) 



(89) 
(90) 



The magnetic field profiles correspond then to a sit- 
uation where the stretching of the magnetic field lines by 
the (imposed) velocity field is compensated by the diffusion 
following equations (|56|l (|57[) . These can be solved if some 
boundary conditions can be applied, which is not obvious 
as the region where the magnetic field is passive does not 
extend to infinity (see Section T4,f .31 for a simple model). Af- 
ter two successive integrations, the radial induction equation 
gives 

fe I .(C) = b r iC--^C 3 , (91) 

where b r i is some unknown constant to be determined by the 
boundary conditions, and we used b r (( — 0) =0 to remove 
one integration constant. Similarly, the azimuthal induction 
equation gives 



MO = vc + - 



brl _ Uj>£\ 

4 



r 3 - v Ur2 c 5 

3 J K a 40 . 



(92) 



where b^i is to be determined by boundary conditions. 



4-1-2 Force-free magnetic field 

When p <C f , nothing can compensate the Lorentz force, so 
the magnetic field has to be force-free: the current is par- 
allel to the magnetic field lines. With our assumption that 
the field is almost vertical, this means that the radial and 
azimuthal currents vanish: 



d(b r 



-D B = 



(93) 
(94) 



Applying the boundary conditions at infinity, we find that 
they are satisfied anywhere in the force-free region: 



b r = ±fe rs + D B (, 



(95) 
(96) 



where ± stands for sgn(£). 

Because the current vanishes, the magnetic field cannot 
diffuse: the velocity is dictated by the fact that the fluid is 
frozen into the magnetic field lines. In particular, isorotation 
is enforced along the magnetic field lines, determining the 
azimuthal velocity to be 



Ut/> 



(97) 



where u'^q is again a constant to be determined by boundary 
conditions. The radial velocity is constant and equals the 
speed u^p at which the magnetic field is being advected or 
diffused. 



In this subsection, we build an approximate model of the 
vertical profiles of velocity and magnetic field by assuming 
that for ( < (b they behave as described in Section 14. f .1 1 
(passive field) and for £ > £g they behave as in Section al .21 
(force- free field). By doing so we neglect the thickness of the 
transition where P is of order unity. To build the model, we 
need to connect the two regions of passive field (£ < (b) and 
force-free field (( > (b), thus determining the proper bound- 
ary conditions at £ = £ b (the conditions at £ = — C,b then 
being given simply by reflectional symmetry). Four condi- 
tions are needed to constrain the four unknowns bri, 6<#>i, 
u'^o and My,. 

Two boundary conditions can be obtained from the 
analysis of the induction equation. The radial component 
states that the azimuthal electric field is independent of the 
height £ (equation I75p . Using equations (|85p and (|9f | at 
( — 0, we obtain a first condition: 



— — (b r l — D B ) + UrO- 



(98) 



The azimuthal component of the induction equation is 
more complicated because b r acts as a source term in this 
equation. By integrating between two heights £i and (2, we 
get the relation 



ra „ , 1 3 



<2 



b r d(. 



(99) 



a 



We can use this relation between and £j to connect the 
passive-field region and the force-free region. As mentioned 
above, we neglect the width of the intermediate region, so 
that the right-hand side of the equation vanishes. Then 



V 



(100) 



where we have used the force-free condition c?c&0(Cb) ~ 0- 
This boundary condition corresponds to a continuous radial 
electric field across C,b- 

Finally, two more conditions are needed in order to con- 
strain the model. These should come from the two compo- 
nents of the equation of motion. Integrating this equation 
across the intermediate region is not straightforward as it 
contains terms proportional to pPo, which vary under our 
simple model from infinity in the passive-field region to zero 
in the force-free region. To build a simple model, we first 
make the simplest assumption that the magnetic field does 
not vary significantly at the transition between the two re- 
gions: 



br(G) = b r ((£) = b rs + D b (e 
b<j>{G) = MCs) = 



(101) 
(102) 



Although not well motivated theoretically, this assumption 
might seem intuitively consistent with the assumption of an 
infinitely thin transition. As shown in the next subsections, 
it is rather well verified for the azimuthal field but not so 
well for the radial field. Indeed, the profile of the radial com- 
ponent of the magnetic field can show a significant variation 
around £b, which is not reproduced by the above descrip- 
tion. The comparison with numerical calculations described 
in the next few subsections, however shows that this model 
is a useful approximate description, which reproduces most 
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relevant effects. A better description of the boundary condi- 
tions is developed in Appendix A for the special case of van- 
ishing viscosity, and heuristically generalised in Appendix B 
to non-vanishing viscosity in order to obtain more precise 
analytical expressions in a two-zone model. 

4.2 Magnetic field diffusion (b rs and Db) 

We recall that our problem is linear and that the full solution 
of the problem consists of a superposition of the response of 
the system to various source terms (see Section 13. 8|) . We 
consider first the source terms b rs and Db, which cause the 
magnetic field to bend and diffuse through the disc. The pro- 
files computed numerically (including only the source terms 
&r-s and Db) are illustrated in Figures U and f5] for three dif- 
ferent strengths of the magnetic field. The quantities that 
are plotted are the dimensionless magnetic field and veloc- 
ity variables as defined in Section 13.11 Note that in each 
case (b rs = 1 and Db = 1) the poloidal magnetic field bends 
outwards. 

The transition between the passive-field region close to 
the midplane and the force-free region at £ > Cs is clearly 
visible for the weak fields (fto = 10 4 and 10 8 ) as the height 
above which the velocity departs from zero and the magnetic 
field tends to the values dictated by the boundary conditions 
at infinity. Note that in the case of fio = 10 the magnetic 
field is sufficiently strong so that it does not really behave 
as a passive field at the midplane. As expected from Sec- 
tion fJTTT] in the passive-field region, the velocity vanishes 
while the radial component of the magnetic field has a linear 
profile and the azimuthal component looks like a third-order 
polynomial. In the force-free region, the azimuthal magnetic 
field vanishes and the radial component is either a constant 
(brs, Figure[4| or linear (DsC> Figure [5]) as required by the 
boundary conditions at infinity. Also as expected, the radial 
velocity is constant and equals the transport velocity of the 
magnetic flux, while the azimuthal velocity is either linearly 
increasing or parabolic. 

Let us now use the two-zone model of Section 14.11 to 
derive approximate vertical profiles. As already mentioned, 
when only these source terms are included, the velocity van- 
ishes in the passive-field region. The radial magnetic field 
profile in this region is then linear and can be obtained us- 
ing equation (1101 [) : 



br = brlC = ( ^ + Dl 
SB 



(103) 



Using equation (|102[) (with b^ — 0), one can obtain the 
azimuthal magnetic field profile arising from the stretching 
of this radial field: 



IV 
4 a 



Cs" 



+ d b ccc 2 - cl)- 



(104) 



The transport velocity of the magnetic flux is then obtained 
from equation (|98|) : 

a brs /,nO 

r C,B 

and the azimuthal velocity profile in the force-free region is 
obtained using equation (|100p : 



u 4> = ( ^ICI - Cs 



+ ^(3 C 2 



At the transition between the two regions, the radial 
and azimuthal components of the velocity show a jump 
qualitatively consistent with the boundary conditions of the 
approximate model (equations 19814100)) . As assumed in the 
model, the azimuthal component of the magnetic field goes 
smoothly to zero at Cs; however, the radial component has 
a behaviour that is not described by the simple boundary 
conditions of Section \4. II Instead of smoothly increasing to 
its final value at infinity, the radial magnetic field reaches a 
maximum and then decreases to its value at infinity. This 
can be qualitatively understood with the following argu- 
ment. At C ~ Cb, the azimuthal velocity has to increase 
and take positive (i.e. super-Keplerian) values, due to the 
requirement that the radial electric current be continuous 
fequation llOOp . This requires a radial Lorentz force directed 
inward, which corresponds to a radial magnetic field decreas- 
ing with £. Given that in the passive-field region the radial 
component of the magnetic field is increasing, it has to show 
a maximum at the transition around £ ~ £b- A quantitative 
description of this effect is given in Appendix A and B. 

To go further in the comparison between the approxi- 
mate analytical model and the numerically computed pro- 
files, we remark that the dependence on the magnetic 
Prandtl number and the magnetic field strength can be 
scaled out of the analytical expressions in the following way 
(here for the source term b rs only): 



br(C<Cs) = 



^7T h& < Cb) 

' SB 

— MCXb) 

a 

7^(C > Cs) 



, c 

brs — , 

Cs 

Ws C (C 



4 Cs VC 



- 1 



— brs 



1 'IS- 1 



(107) 
(108) 
(109) 
(110) 



CI) 



(106) 



In Figures [6] and [7] we show the scaled profiles (i.e. of b r , 
ab r f > /(V(, B ), VC,BU r /a, and w^/Cs) as a function of C/Cs for 
different values of the magnetic field strength (Figure[S} and 
of the magnetic Prandtl number (Figure [7J|. The model sug- 
gests that all these curves should lie on top of each other, 
and coincide with the dashed curves that represent the ana- 
lytical profiles given by equations (|107|) - (|110l) . The different 
profiles do indeed lie very close to each other (except maybe 
for the /3o = 100 case, where the passive-field hypothesis 
at the midplane is less well verified) showing that the de- 
pendences on the magnetic field strength and the magnetic 
Prandtl number have been properly scaled out. 

One exception to this agreement, however, is the bump 
in radial velocity close to C — Cs, which is more and more 
pronounced as /3o increases. The two-zone model fails to 
reproduce this feature as it takes place in the intermedi- 
ate region where the magnetic field is neither negligible nor 
dominant. Furthermore, note that the agreement between 
the numerical profiles and the analytical ones is not perfect. 
All the differences can actually be traced back to the profile 
of radial magnetic field, which goes through a maximum at 
C ~ Cs as discussed earlier. The presence of the maximum 
leads to a larger radial field at C = Cb th an assumed in the 
model. The bump in the radial velocity is due to the pres- 
ence of the maximum in the radial magnetic field, together 
with the requirement that the azimuthal electric field is uni- 
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Figure 4. Vertical profiles of the radial (left) and azimuthal (right) magnetic field (top) and velocity (bottom) in response to a radial 
component of the surface magnetic field b rB = 1 (i.e. the solution vector Xt, rs ). The different colours correspond to different strengths 
of the vertical magnetic field: /3o = 10 (red), /3q = 10 4 (blue), and /3q = 10 s (black). The vertical dotted lines illustrate the height £g at 
which the magnetic and thermal pressures coincide, marking the transition between the passive and the force-free field regimes. 




Figure 5. Same as Figure [4] but for the source term Db = 1 describing the radial gradient of the vertical magnetic field strength (i.e. 
the solution vector Xpg). The dashed and dot-dashed lines show the prediction of the analytical model respectively with the simple 
boundary conditions of Section l4.1l and with the more precise boundary conditions of Appendix B (only for the case /3q = 10 8 ). 



form. The improved boundary conditions described in the 
Appendix B yield a much better agreement between the nu- 
merical profiles and analytical ones (shown with dot-dashed 
lines in Figures [6] and O . 

Figure [8] shows the dependence of the transport velocity 
of the magnetic flux on the strength of the vertical field 
(upper left and right panels for b rs and Db respectively). 



The prediction of the analytical model agrees well with the 
numerical calculation for /?o > 10 3 — 10 4 , especially when 
using the improved boundary conditions of Appendix B. We 
recall that for a thin disc and a significant inclination of 
the surface magnetic field, the dominant contribution to the 
diffusion is that of the radial field ut rB (which is larger than 
udb by a factor ~ r/H). It is interesting to note that this 
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Figure 6. Same as Figure [4] but rescaled according to the analytical model of Section |4. II The different colours correspond to different 
strength of the vertical magnetic field : /3q = 10 2 (red), /3rj = 10 4 (blue), /3o = 10 s (black). The dashed and dot-dashed lines show the 
prediction of the analytical model respectively with the simple boundary conditions of Section 14.11 and with the more precise boundary 
conditions of Appendix B. 




diffusion process is less efficient for weak magnetic fields. 4.3 Transport due to angular momentum loss in 

Indeed for the parameters explored here, the diffusion is up an outflow (6</> s ) 

to a factor of 4 slower than the crude estimate stated in the 
introduction (equation J2J, which corresponds to ut rs = 1 in 
Figure [8]). This slower diffusion is a consequence of the fact 

that field lines bend on a length-scale Cs, which is larger for Applying the two-zone model with the simple boundary con- 
weak fields. ditions of Section 14.11 predicts that no advection is induced 

by the source term b<p B . Indeed, in the passive-field region the 
velocity vanishes, as does the radial magnetic field. Alone, 
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Figure 8. Contributions to the transport velocity of the magnetic flux, u^, as a function of the magnetic field strength, for different 
source terms: diffusion due to the bending across the disc (b rB source term, upper left panel), diffusion due to the radial gradient of B z 
(Db, upper right panel), advection due to the turbulent viscosity ('hydrodynamic' source terms, Dh = 1 and D„e = 0, lower left panel), 
advection due to angular momentum removal by an outflow (6^, a , lower right panel). In all panels, the full black line corresponds to the 
numerical results, while dashed and dash-dotted lines correspond to the analytical model respectively with the simple boundary conditions 
of Scction l4,f l and with the more precise boundary conditions of Appendix B. The red lines show the corresponding contributions to the 
mass transport velocity u m - 



the azimuthal magnetic field has a linear dependence on (: 



Cb 



(111) 



The transport of magnetic flux vanishes: = 0, while the 
azimuthal velocity jumps around C,b to its final value: 



(112) 



This is not the whole story however, because two effects 
have not been taken into account. First, the jump in radial 
magnetic field around £s (Figure [9| induces some trans- 
port of magnetic flux. This is described quantitatively in 
Appendix B, and the analytical prediction for the transport 
velocity of magnetic flux is compared to numerical results in 
Figure [8] The agreement is very good for rather weak fields 
(/3o > 10 2 ). For larger values of the magnetic field, another 
effect neglected in the two-zone model becomes important. 
Indeed, the magnetic force cannot be neglected anymore and 
the torque exerted by the azimuthal field induces a non- 
negligible radial velocity, which contributes to the transport 
of magnetic flux. For fio < 100, the transport velocity of 



magnetic flux is then very close to that of mass (given by 
equation (|71[> ). 

As might have been expected, the efficiency of the sur- 
face azimuthal magnetic field to transport mass and mag- 
netic flux strongly increases as the magnetic field strength 
is increased (Figure [8}. It is interesting to note however that 
the transport velocity of magnetic flux is much larger than 
that of mass, when the magnetic field is rather weak. 



4.4 Advection due to the turbulent viscosity 

All the remaining source terms {D v t,, Dh and those which 
solution vector is Xk) are hydrodynamical in nature and 
describe the advection process due the turbulent viscosity. 
First, we note that the source term D„s admits an exact 
solution with uniform radial and azimuthal velocities and 
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Figure 9. Same as Figure [4] but for an azimuthal surface magnetic field = 1. Note that, in reality, b^s would be expected to be 
negative (so that the outflow removes angular momentum from the disc), so this component of the solution would be multiplied by a 
negative coefficient. 




vanishing radial and azimuthal magnetic fields: 

b r = 0, 
b<t, = 0, 
u r = — 3aD u s, 

2 



(113) 
(114) 
(115) 

(116) 



We now focus on the other hydrodynamical source 
terms, and their solution in the approximate analytical 



model. In the passive-field region the velocity profiles are 
the parabolae given by equations (1851) (190 1) (note that the 
parabolic shape is clear in Figure 1 10 p . The magnetic field 
profiles result from the stretching of the magnetic field lines 
by the parabolic velocity profile, and from the condition that 
both radial and azimuthal components vanish at £ = £b (for 
the simple boundary conditions of Section [4. ip : 



W = ^-C(Cs - C 

a 6 



(117) 
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h = — (Cb 

a 



C 2 ) 



(118) 



Then the velocities in the force-free region and, most impor- 
tantly, the advection speed of the magnetic flux are found 
using equations and (|100l) : 



. U r2 2 
UrO + ~Cs, 



(119) 



(120) 



%>(C > Cs) = ""00 + ^-Cs + -TT—(b- 

3 15 a 

These profiles are compared with the numerical solution in 
Figure [10] for the case /3o = 10 8 (the analytical profile is the 
dashed line). The agreement is acceptable though the mag- 
netic profiles are slightly underestimated. Again, a quanti- 
tatively better agreement is found when using the improved 
boundary conditions of Appendix B (dot-dashed lines). 

Rewriting the transport velocity of the magnetic flux as 
a function of the source terms, one finds 



- ~ + 5D H -3A,s + ^~ 
a 2 4a 2 



U,<j, 



5D t 



r 2 

2 , C,B 



4a + 



(121) 



This shows that the inward transport is faster for weak fields 
(i.e. large £b) if Dh > 0.6. It is for example the case for our 
fiducial values of Dh = 1 and D„£ = 0: 



4q,2 



r 2 

Aa 2 + ^ 



(122) 



The lower left panel of Figure [8] compares the transport ve- 
locity of the magnetic flux Mhyd (full black line) with equa- 
tion (I122p (dashed line) and the transport velocity of mass 
(red line). Interestingly, the transport velocities of mass and 
magnetic flux coincide for strong fields /3o < 10, but for weak 
fields the advection of magnetic flux is much faster than that 
of mass (by a factor up to 10 for the present parameters). 
This behaviour is correctly reproduced by the approximate 
analytical model, and its interpretation will be discussed in 
the following subsection. 



4.5 Advection/diffusion speed as a vertical 
average 

The results of Sections l4.2l l4.3l and l4.4l show that the trans- 
port of magnetic flux is significantly affected by the vertical 
structure, in the sense that it differs from the s imple vertical 
averag e used by many previous works such as iLubow et al.l 
(1994). This can be understood by noting that the trans- 
port velocity of magnetic flux, contrary to that of mass, is 
not given by a density weighted average. B elow, we repeat 



with our notation the argument given by lOgilvie fc Livid 
|200lf ) on the correct way of averaging the equations in or- 
der to obtain the effective transport velocities of mass and 
magnetic flux. We then use it to interpret our results. 

Dividing equation (|75[) by the dimensionless resistivity 
a/V, and integrating between £ = and £, we obtain: 



(123) 



where r\ is an average dimensionless resistivity defined in 



the following way (it is actually the inverse of the height- 
averaged conductivity) : 

fj = — ^ , (124) 

and u r is the average of the radial velocity weighted by the 
conductivity (1/ 77) : 

Ur= jh^l!^ dC (125) 

This averaging procedure shows that the contribution of 
advection to the velocity at which the magnetic flux is trans- 
ported is the vertically averaged radial velocity weighted by 
the electrical conductivity. This average may be very dif- 
ferent from the density-weighted average u m that describes 
mass transport. 

Note that equation (1123[) is valid for any height £ up 
to which the averaging procedure is performed, and which 
can be chosen in the most convenient way. Integrating up 
to an infinite height may not be very useful as the force-free 
regime would dominate the average and one would simply 
recover = Ci, as was already found in Section [2] Instead, 
the two-zone model developed in Section [4.11 suggests that 
performing the average up to the height (b is most useful. 
By doing so, one can indeed recover the results for u^> of 
Sections 14.21 and 14.41 in a simpler way: 



-D, 



+ u r 



(126) 



where u r is the average advection velocity as defined by 
equation (|125[) . performing the averaging procedure between 
C, — and ( = ( B , and we used f\ — a/V since the diffusion 
coefficients are uniform. Following the two-zone modelling of 
Section l4.1l the average advection velocity can be computed 
using the purely hydrodynamical velocity profile (because 
the field is assumed to be passive for £ < £ B ): 
Cb 

a, - / (UrO + u r 20 d( = UrO + : ^Cr- (127) 



/ a2\ j j- , U r2 t 1 
(UrO + U T 2Q ) AC, = U r H ^~Cs- 



1 

Cb j 

u r gives the value of the advection speed due to the hydro- 
dynamic source terms, as described by equation (|165p . The 
value of b r at ( B then has to be prescribed by the bound- 
ary conditions. Using that of Section 14.11 gives ^(Cb) — 
b rs + DbC,b, and one recovers equation (JTUSJ) . Using the 
boundary conditions of Appendix B is also possible in prin- 
ciple but it requires first solving for the full vertical profiles. 

The above argument allows for a deeper interpretation 
of the higher advection speed of weak fields, which was ob- 
served in the previous section. Indeed, the mass transport 
velocity u m corresponds to a density-weighted average and is 
therefore dominated by the velocity close to the midplane, 
where the density is large. By contrast, the magnetic flux 
transport velocity u^ is a conductivity-weighted average, 
which for the case of uniform diffusion coefficients is equiv- 
alent to a height-average. As a consequence the advection 
of a weak field is significantly affected by the low-density 
region away from the midplane. This low-density region at 
high £ can have a much faster inward radial velocity than in 
the midplane, if it r 2 < as is readily found for standard pa- 
rameters. Such parabolic velocity profiles, where the radial 
velocity becomes more and more negative with the distance 
from the midplane, are well kn own solutions of visco us hy- 
drodynamical disc models (e.g. iTakeuchi fc Linll2002l ). 
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Figure 11. Profiles of the radial (left) and azimuthal (right) components of the magnetic field (upper) and the velocity (lower), in a 
stationary situation where the advection of magnetic flux is exactly compensated by the diffusion of the radial field across the disc. The 
parameters used are: 0g = 10 4 , 'hydrodynamic' source terms Djj = 1 and = 0, and a radial surface magnetic field b rB ~ 14.5 chosen 
so that diffusion compensates advection (note, however, that our model neglects the magnetic compression of the disc). 



4.6 Inclination of the field lines 

As outlined in the introduction, one of the most important 
diagnostics of the relative efficiency of the magnetic flux 
advection and diffusion processes is the inclination of the 
field lines at the surface of the disc that can be achieved 
in a stationary situation. In such a stationary situation the 
diffusion velocity due to the inclination, UbrsKs (recall that 
b rs = r/H x Brs/Bz), compensates the advection velocity 
ithyd (as an example, let us focus on the case Dh = 1, D„s = 
0). Figure[TT]illustrates the vertical profiles obtained in such 
a stationary situation. The inclination is then: 

B^ = _Hu^ 
B z r u brs 

In the limit of an infinitely thin disc, and if Uhyd and ut rs are 
of order unity, one recovers the classical result that the incli- 
nation is vanishingly small. The inclination of the field lines 
can be estimated with the approximate analytical model 
(here with the simple continuous magnetic field boundary 
condition of Section 14. ip : 



B rs _ H 
~TT~ - —PQb 

B z r 



Aa 2 



4<r 



t 2 



(129) 



In Figure [121 this estimate (represented by dashed lines) is 
compared to the numerical result. The more precise estimate 
represented in dash-dotted line is given in Appendix B. As 
expected, the formula shows that the bending does not de- 
pend much on the value of a (at least if a 2 <C 1) but rather 
on the magnetic Prandtl number V ■ In particular one recov- 
ers the result that a significant inclination can be achieved 
if "P ~ r/H (but this is not expected to be true). The for- 
mula and associated figure however suggest a new way to 
achieve a large inclination: to be in the weak field regime 
where /3q and £g are large. Indeed, a significant inclination 



is possible for weak fields if the disc is not too thin. For ex- 
ample, if H/r = 0.1 a field of /3o ~ 10 3 can be significantly 
bent (B r ~ B z ). This is due both to a reduced diffusion 
and a faster advection speed as compared to the rough esti- 
mates quoted in the introduction (as discussed in previous 
subsections). 

Similarly the ratio of the radial and azimuthal compo- 
nents of the surface magnetic field can be estimated in a 
stationary situation where diffusion is compensated by the 
advection due to an outflow: 

B r s = _W6^ (13fj) 

J3<j)s Ubrs 

(Recall that B^ s is negative so that the outflow removes an- 
gular momentum from the disc) . This ratio is represented in 
the lower panel of Figure[l2]and compared with the estimate 
of the analytical model given in Appendix B. As might be 
expected, the advection due the outflow is more efficient at 
bending the magnetic field for strong magnetic fields. 

In order to relate the above ratio to a radial inclination 
of the field with respect to the vertical, one would need to 
know the value of the azimuthal magnetic field at the sur- 
face, which determines the efficiency with which the outflow 
removes angular momentum. This would require a study of 
the launching and acceleration of the outflow, which is be- 
yond the scope of this paper. We therefore keep it as a free 
parameter and obtain the following expression for the incli- 
nation: 

B rs _ Bc/ys U60 s 
B z B z Ubrs 

For a given value of B^/Bz, Figure [P2l suggests that only 
magnetic fields above a certain strength can sustain a sig- 
nificant inclination of the field lines, if the advection is due 
to an outflow. For B^s/B z = —1, the threshold is /3o ~ 10, 



(131) 



Magnetic flux transport in accretion discs 19 



80 



60 



40 



20 









/ 






/ 


/ ' 
/ J 

/ A 
// 






/ 

/ y 
' //' 





























10 u 



Oa. 1 



10 U L 
0.01 



H/r 



I 0.00 



1 .00 r ' 



0.10 



0.0 




10 4 



10 b 



10 b - 



10 4 



1 o" - 



I 0.0 



I 00.0 



Figure 12. Radial inclination of the field lines at the surface 
of the disc obtained in a stationary situation. The upper panel 
shows the inclination (B rB /B z ) in units of H/r such that the 
diffusion compensates the advection due to the turbulent viscosity 
('hydrodynamic' source terms). The lower panel shows the ratio 
—brs/bfa such that diffusion compensates advection due to an 
outflow, fn both panels, the full black line corresponds to the 
numerical results, while dashed and dash-dotted lines correspond 
to the analytical model respectively with the simple continuous 
magnetic field boundary conditions (Section 14.11) and with the 
more precise boundary conditions of Appendix B. Note that in 
the lower panel the dashed line is not shown as it would equal 0. 



Figure 13. Upper panel: maximum strength of the magnetic 
field (i.e. minimum 0q) for which an inclination of 30° of the 
surface magnetic field can be sustained, if the advection is due 
to the turbulent viscosity. This maximum strength is plotted as 
a function of the aspect ratio of the disc H/r. Bottom panel: 
minimum strength of the magnetic field (i.e. maximum /?rj) for 
which an inclination of 30° of the surface magnetic field can be 
sustained, if the magnetic flux advection is due to an outflow. 
This minimum strength is shown as a function of the azimuthal 
magnetic field component at the surface of the disc, parametrising 
the efficiency at which the outflow extracts angular momentum 
from the disc. 



while for B^ s /B z = -10 it reaches /3 Q ~ 10 3 . Although 
magneto-centrifugal acceleration is usually thought to pro- 
duce an outflow if the magnetic field is strong (/Jo of order 
or slightly above unity), we note that outflows with weaker 
fields may be driven by discs wi th superheated surface lay- 
ers. Indeed. [Murphy et alj l|2010t) have recently reported the 
production of a magneto-centrifugal outflow from a weakly 
magnetised disc (/Jo ~ 100 — 1000). This opens the possi- 
bility that an outflow from a weakly magnetised disc could 
significantly affect the transport of magnetic flux. 

4.7 A scenario to obtain strong inclined magnetic 
fields 

The above discussion suggests a scenario by which strong 
and inclined magnetic fields could be obtained in a moder- 
ately thin disc through the inward advection of an initially 
weak magnetic field. In such a disc, the advection of a weak 



field by the accretion flow due to the turbulent viscosity is 
possible up to a certain strength, which depends strongly 
on the aspect ratio of the disc (Figure 1131 upper panel). 
Unless the disc is extremely thick, this strength is however 
rather weak: for example /Jo ~ 10 3 is reached for an already 
rather thick disc with H/r > 0.07. Creating and sustaining a 
stronger inclined field might then be possible if the advection 
of the magnetic flux is primarily due to an outflowQ, which 
is possible above a certain strength of the magnetic field. 
The minimum strength allowing the outflow to compensate 
diffusion is shown as a function of B$ s /B z (parameterising 
the ability of an outflow to remove angular momentum from 
the disc) in Figure [13] (lower panel) . 

2 Note that this does not mean that the angular momentum 
transport and hence the advection of mass is primarily due to 
the outflow, since the advection speed of mass and magnetic flux 
can be very different. 



20 Jerome Gullet and Gordon I. Ogilvie 



1 .00 




-B,/B z 

Figure 14. Aspect ratio above which a disc can accumulate mag- 
netic flux and create strong magnetic fields in the scenario out- 
lined in Section 14.71 This aspect ratio depends strongly on the 
efficiency at which an outflow removes angular momentum, pa- 
rameterised by —B^ s /B z . 

A condition for this scenario to work is that the max- 
imum strength for advection through turbulence should be 
larger than the minimum strength for advection through an 
outflow. This requirement is fulfilled above a certain aspect 
ratio, shown as a function of B^/Bz in Figure 1141 This 
suggests that for reasonable values of the surface azimuthal 
field, the disc needs to be rather thick to enable the produc- 
tion of a strong magnetic field: H/r > 0.2 for B# s /B z = — 1, 
or H/r > 0.07 for B^/B z = -10. 

Obviously, these numerical values should be taken with 
a grain of salt for several reasons. The effect of an outflow 
has been taken into account in a quite approximate way. A 
proper description of this scenario would need to include the 
launching and acceleration process of the jet, which should 
determine self consistently the value of the azimuthal com- 
ponent of the magnetic field at the surface of the disc. Fur- 
thermore a variation of the diffusion coefficients in the ver- 
tical direction might be expected to change quantitatively 
the result, as will be studied in paper II. The turbulence 
may also not have such a simple effect as isotropic diffusion 
coefficients ; for example, the vertical diffusion may be ex- 
pected to be less efficient than the radial one due to smaller 
vertical velocities in MRI turbulence. As a consequence, the 
numerical values for the minimum aspect ratio should not be 
considered as a strong constraint, and it might turn out that 
this scenario is effective for somewhat lower aspect ratio. 



5 CONCLUSION AND DISCUSSION 

We have developed a formalism enabling us to compute self- 
consistently the radial transport of poloidal magnetic flux 
and mass in a thin accretion disc, where the effect of tur- 
bulence is modelled by an effective viscosity and resistiv- 
ity. For this purpose, we performed a systematic asymp- 
totic expansion in the limit of a thin disc. We assumed that 
the magnetic field is only weakly inclined with respect to 
the vertical direction, in order to avoid the complication of 
modelling a magneto-centrifugal outflow, as well as for self- 
consistency in the expansion. Owing to this assumption, the 



solution (and thus the transport velocity of the magnetic 
flux) depends linearly on a number of source terms. These 
source terms are the radial gradients of various quantities 
(surface density, aspect ratio, magnetic field strength), and 
non-vanishing values of the radial and azimuthal magnetic 
field at the surface of the disc. Each of them describes a dif- 
ferent physical process which can then be studied and under- 
stood independently from the others: the diffusion due the 
bending of the field lines across the disc (b ra ), the diffusion 
due the radial gradient of its intensity (Db), the advection 
due to angular momentum loss in an outflow (6^ a ), and that 
due the turbulent viscosity ('hydrodynamic' source terms). 

An important aspect of this calculation is that it takes 
into account the vertical structure of the disc. The vertical 
profiles of velocity and magnetic field are indeed explicitly 
computed. The formalism also allows for a vertical depen- 
dence of the diffusion coefficients. In this paper, we have 
however only considered the idealised case of uniform coef- 
ficients, leaving the question of their vertical structure for 
paper II. Even in this simplest case, we found the vertical 
structure to have an important impact on the transport of 
magnetic flux. Indeed, the simple vertical averaging usually 
performed gives results that can be in error by a factor of 
10 or more, because of its improper treatment of the verti- 
cal structure. We have developed an approximate analytical 
model which is able to describe this behaviour, and provides 
a physical explanation. 

Two main effects of the vertical structure on the trans- 
port of magnetic flux have been uncovered by this study. 
Both of these effects are important when the poloidal mag- 
netic field is rather weak in the sense that the magnetic 
pressure is small compared to the thermal pressure in the 
midplane. First, we showed that the diffusion of weak fields 
is less efficient because the bending of the field lines takes 
place on a larger vertical scale. Second, the advection of 
magnetic flux can be significantly faster than that of mass, 
owing to fast radial velocities in the low-density regions away 
from the midplane. 

Taken together, these two effects enable a more efficient 
inward transport of magnetic flux, which partly alleviates 
the long-standing problem of too fast outward diffusion of an 
inclined magnetic field. Indeed, they suggest that in a mod- 
erately thin disc, magnetic fields below a certain strength 
can be advected inwards without significant outward diffu- 
sion. This critical strength depends steeply on the aspect 
ratio of the disc, however, and is probably rather weak for 
reasonably thin discs. We propose a scenario in which the 
advection due to the presence of an outflow could then do the 
rest of the job, and provide a means to obtain yet stronger 
magnetic fields. A condition for a successful advection is that 
the aspect ratio of the disc be larger than a critical value, 
which is dependent on the efficiency of the outflow to extract 
angular momentum from the disc. Confirming or invalidat- 
ing this scenario will require a better understanding of the 
magnetic flux advection induced by an outflow with a weak 
magnetic field. 

The faster advection of weak magnetic fields could also 
have interesting consequences on the time-dependence of 
magnetised accretion discs. Indeed, the magnetic flux dis- 
tribution would have the counterintuitive ability to evolve 
significantly faster than the mass distribution. These two 
potential consequences should be studied in a global time- 
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dependent disc model, where the inclination at the disc sur- 
face is computed self-consistently as a result of the distribu- 
tion of magnetic flux. The time-evolution of the mass and 
magnetic flux could then be described using the transport 
rates computed with the formalism described in this paper 
(with the equations given in Section [2]) • This question is left 
for future work. 

One limitation of our analysis is that it is restricted to 
small values of the inclination of the magnetic field with 
respect to the vertical direction. By comparing with di- 
rect numerical simulations of a stratified shearing box, we 
found that the small-inclination approximation gives good 
results even at significant inclination as long as an outflow 
is not launched (i < 30°) and the field is not very strong 
(/?o ~ 10 2 — 10 4 ). However, it is expected to break down at 
stronger fields because magnetic compression would become 
important, or at still larger inclination, in which case the 
launching process of an outflow should be included. 

A related worry could come from the large values of 
obtained in some cases for weak magnetic fields (e.g. Fig- 
ure [TO]). If the aspect ratio is only moderately small (e.g. 
H/r ~ I0~ 2 — 10 _1 ) then the azimuthal component of the 
magnetic field may be larger than the vertical one. In ad- 
dition to contradicting the initial assumption (in particular 
the neglect of compression), one might worry about the sta- 
bility of such a solution. 

We have underlined in this paper the importance of the 
vertical profile of radial velocity. iFromang et al.l l|201ll ) have 
studied this profile by means of global MHD simulations and 
found that it differed from the prediction of viscous alpha 
models. This comes from the different vertical profile of the 
turbulent stress, compared to a viscous model with a viscos- 
ity independent of height (as was assumed in this paper). 
This shows the limitation of the idealised model studied in 
this paper. In paper II, we will consider a vertical profile 
of the viscosity in order to obtain a more realistic vertical 
profile of the stress. This study confirms the robustness of 
our main result: that the magnetic flux is advected faster 
than mass when the field is weak. 

Both the faster advection and the reduced diffusion of 
weak magnetic fields rely on the dynamics of the low-density 
region away from the midplane. Turbulent MHD simula- 
tions both of stratified shearing boxes and global mod- 
els show that this region tends to be strongly magnetised 
jMiller fc Stondl200d ; lFromang fc Nelsonlbood i. and it is of- 
ten called a 'corona'. (Note that 'strongly magnetised' here 
refers to the fluctuating component of the magnetic field 
which has a pressure comparable to the thermal pressure, 
but this does not mean that the large-scale poloidal mag- 
netic field is strong.) This article thus suggests that the dy- 
namics of the corona and its connection with the transport of 
magnetic flux should be studied in more detail. It is unclear 
whether this corona should be modelled with an effective 
viscosity and resistivity like the main part of the disc, since 
some authors refer to this corona as non-turbulent due to its 
strong magnetisation. The vertical extension of the corona 
would also be of some interest for the diffusion process. In- 
deed, a more extended corona (because it is hot or because 
of magnetic support) could help reduce even more the dif- 
fusion of weak fields by increasing the length-scale on which 
the field lines are bent. 

Finally we should stress that all this analysis relies on 



the assumption that the turbulence can be modelled by an 
effective resistivity and viscosity. It is not clear that this 
assumption is valid, especially for diffusive processes acting 
on a scale comparable to the scale-height. Indeed, in that 
case there might not be a real scale separation since the 
correlation length of the turbulence could be of the same 
order of magnitude. In principle then, the right method 
would be to perform turbulent numerical simulations that 
self-consistently describe the effect of turbulence. MHD sim- 
ulations of a stratified disc containing a significant net ver- 
tical flux have however proved problematic so far. To avoid 
the further complication and the necessarily low resolution 
of global simulations, it would be useful to study the ad- 
vection and diffusion of the magnetic flux in a local model. 
The classical stratified shearing box model would need some 
modification in order to take into account the relevant source 
terms. The inclusion of the non- vanishing horizontal compo- 
nent of the magnetic field at the surface of the disc could be 
achieved simply by changing the vertical boundary condi- 
tions, as has been done in this article for ID non-turbulent 
simulations. A boundary condition that imposes the mag- 
netic field may be problematic in the presence of turbulence 
however. The inclusion of other source terms describing the 
advection due to the effect of turbulence is much less obvi- 
ous, but may deserve to be further investigated. 
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APPENDIX A : A DESCRIPTION OF THE 
INTERMEDIATE REGION IN THE CASE OF 
ZERO VISCOSITY 

In this appendix we provide an analytical description of 
the intermediate region around £b connecting the passive 
field and the force free field regimes. For this description 
to be tractable we assume this intermediate region to be 
thin, and we consider only the case of vanishing magnetic 
Prandtl number (i.e. zero viscosity a, but finite resistivity 
Q./'P). We concentrate on the source terms Db and b rs in 
order to explain the bump observed in the profile of the 
radial component of the magnetic field (the hydrodynami- 
cal source terms with vanishing viscosity would not be very 
meaningful anyway). 

The density profile around C,b is approximated by an 
exponential: 



where we defined x = £ — (b- This expansion of the density 
profile is valid if x <C C,b ■ Note that the intermediate region 
where /3 ~ 1 has a typical thickness of x ~ 1/Cbi therefore 
the assumption of a thin transition is valid if 3> 1- The 
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differential system is then rewritten in the following way: 



U4 , + e <BX (d x b r - D B ) 



d x ^^(d x b r ~ D B ) + u r 
d x \^d x b$ 



u 4>\ = 



0, 
0, 

0, 

-b r . 



(133) 
(134) 

(135) 
(136) 



The term on the right hand side of Equation (|136p is further- 
more neglected due to our assumption that the transition 
region is thin. Equations (I135p - (|136p then simply state that 
the azimuthal and radial electric field are constant, and their 
value can be estimated at (bi — > ±oo (i.e. at £ = Cb)' this 
is consistent with the boundary conditions given by equa- 
tions jnij and (flUD)) . Thus 



V 



(dxb r - D B ) 



Uij, 



Au<f, — —dxbj,, 



(137) 
(138) 



where = Ur(Cs) ~ "^/^t^cMCe) ~ Db) and Au^ = 
u 4>(Cb) ~ a /^ , ^C^0(Cs)- Introducing these relations into 
equations (|133[) - (|134[) . one obtains the system of two equa- 
tions: 



V 



dxbd 



(dxbr 



V 



(d x b r - D B ) + 4e 



Db) = Am^, 



CBX dxb4, = 



(139) 
(140) 



which can be solved to obtain explicit functions for d x b r and 
d x b^: 



dxb r = D B + ^ 



My 



■ 4e <BX A 



W0 



+ 4e 2 ^ 3 



d x bd, 



£ + 4e 2 ^* 



(141) 
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Integrating over x, one can compute the jump in b r across 
the intermediate region: 



b r (x)~ b r (— x) 



Ttt AMj 
a ( B 



C, b x — In 



2V 



V 

a ( B 



+2D B x, 
(143) 



if e^ BX 3> 1. The two terms proportional to x represent 
the smooth variation in the passive field regime (from —x 
to 0) and force free field regime (from to i), and do not 
correspond to a jump in the transition region, thus 



Ab r = b r (^)-br(Q) = -^^-ln 

a Q B 

Similarly, one obtains the jump in b^ 



2V\ Vu^ 
a J a ( B ' 



Abj, = b^Cs) ~ MCb) 



Vtv U,p 
4a ( B 
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These new boundary conditions can then be injected 
into the two-zone model. It turns out that, in the limit £fj 3> 
1, only one term is significant, and the jump conditions can 
be simplified to 



Ab r 
Ab$ 



Vn Au^ 
a (b 



(146) 
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This comes from the fact that b r /b^, — 0(1/ C/ B )- The two- 
zone model then gives the following slope of the profile of 
the radial component of the magnetic field: 



b r i 



2a v J 



+ Di 



(148) 



Using the same two-zone description for the marginal 
antisymmetric MRI mode, one can compute the value of the 
magnetic diffusivity under the marginal stability hypothe- 
sis. For this purpose, one assumes the mode to be antisym- 
metric about the midplane (contrary to the calculation in 
Section m>- We obtain 



(149) 



a 3tt 
V = ~2~' 

which is in good agreement with the numerical calculations 
at large /3o and without viscosity (respectively dashed and 
dotted line in Figure [2}. Using this value of a/V in equa- 
tion (|148p , one finally obtains the slope of the radial mag- 
netic field profile: 



Wl = 



3 f brs 



+ Di 



(150) 



which reproduces the larger slope observed in numerical pro- 
files in Section 14.21 Note that equation (|148|l indicates that 
relaxing the marginal stability hypothesis would change the 
slope of the radial magnetic field profile to a different factor 
than 3/2. It even shows that if a/V < it/2 the slope and 
hence the transport rate would change sign ! This is not 
believed to be physical, and is probably an artefact of con- 
sidering a stationary flow that is unstable to the MRI. This 
assertion is supported by the fact that the value a/V = tt/2 
corresponds to the largest-scale MRI mode symmetric about 
the midplane being marginally stable (this can be obtained 
using the two-zone model). At still lower values of a, the 
numerical profiles show multiple bending of the field lines, 
which is also an arte fact due to the fact that the flow is 
unstable to the MRI (|Qgilvie fe Livioll200ll '). 



APPENDIX B : TWO-ZONE MODEL WITH 
MORE PRECISE BOUNDARY CONDITIONS 

Here we assume that the boundary conditions found in 
Appendix A can be generalised to non-vanishing magnetic 
Prandtl number in the following way: 

Abr = (151) 
O C,B 

A6 = 0. (152) 



These boundary conditions are the same as equations (|146p 
and ()147p with the value of a/V given by the marginal sta- 
bility analysis ( equation I149[) . This heuristic generalisation 
is motivated by the numerical results, which show that the 
jump in radial magnetic field is independent of the magnetic 
Prandtl number (see Figure 0. This feature is well repro- 
duced by the above boundary conditions. Another way of 
looking at this generalisation is to assume a more general 
form inspired by equations (|146p and (| 147p : 



A6 r 
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a (b 
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where A is some unspecified numerical factor which may 
depend on V. Using then a two-zone model to describe 
the marginal MRI mode leads to the following value: A — 
2a/ (3V). Introducing this value into equation 1|153[) gives 
equation (|151[) . 

We now use the new boundary conditions in conjunction 
with the two-zone model of Section ^. ll to construct analyti- 
cal profiles corresponding to the different source terms. The 
analytical profiles for the source terms b rs and Db are then 
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2 VCs 
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Those for the source term 6^, s are 
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2V Cfl 

And finally the analytical profiles for the hydrodynamic 
source terms are the following: 

2 
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From these, one can deduce the inclination of the field 
lines at the surface of the disc, when the diffusion is com- 
pensated by advection due to effective viscosity (here for 
D H = 1, and D uY . = 0): 
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as well as the ratio of the radial to azimuthal components of 
the surface magnetic field, when the diffusion is compensated 
by the transport due to an outflow: 

Brs Ub6s 2 a 1 
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